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ABSTRACT 

We present a new method for analyzing ion, or 
molecule, distributions around helical nucleic acids 
and illustrate the approach by analyzing data derived 
from molecular dynamics simulations. The analysis 
is based on the use of curvilinear helicoidal coor- 
dinates and leads to highly localized ion densities 
compared to those obtained by simply superposing 
molecular dynamics snapshots in Cartesian space. 
The results identify highly populated and sequence- 
dependent regions where ions strongly interact with 
the nucleic and are coupled to its conformational 
fluctuations. The data from this approach is pre- 
sented as ion populations or ion densities (in units of 
molarity) and can be analyzed in radial, angular and 
longitudinal coordinates using 1D or 2D graphics. It 
is also possible to regenerate 3D densities in Carte- 
sian space. This approach makes it easy to under- 
stand and compare ion distributions and also allows 
the calculation of average ion populations in any de- 
sired zone surrounding a nucleic acid without requir- 
ing references to its constituent atoms. The method 
is illustrated using microsecond molecular dynam- 
ics simulations for two different DNA oligomers in 
the presence of 0.15 M potassium chloride. We dis- 
cuss the results in terms of convergence, sequence- 
specific ion binding and coupling with DNA confor- 
mation. 

INTRODUCTION 

How ions interact with DNA and what role they may play in 
modulating the structure and interactions of the DNA dou- 
ble helix has been the subject of many experimental and the- 
oretical studies recent years. In structural terms, monova- 
lent ions have been the subject of controversy because they 
are difficult to distinguish from water molecules in crystallo- 
graphic studies (1-4) (even at very high resolution, (5)) and 
they are also not amenable to nuclear magnetic resonance 
(NMR) investigations. Thus, how frequently ions bind to 



specific DNA sites is still open to question, although re- 
sults using other ions, which are less biologically relevant, 
but easier to locate experimentally (such as thallium or am- 
monium which are both reasonable models for potassium 
cations) have shown that monovalent ions can bind both in 
the major groove (favoring GC-rich regions) and in narrow 
minor grooves (e.g. A-tracts) (6,7). Experimental efforts in 
this field continue to develop and the reader is referred to 
recent work using small-angle X-ray scattering (8) and so- 
called 'ion counting' spectroscopic approaches (9,10). 

In view of the experimental difficulties, molecular dynam- 
ics (MD) simulations have been used to study ion binding 
for many years (1 1-18). Initially, such simulations were lim- 
ited to a few nanoseconds. Ion penetration of the grooves, 
and, in particular, replacement of waters forming the minor 
groove spine in A-tracts, was observed and was instrumen- 
tal in encouraging further studies, however the relatively 
slow diffusion of the ions made it difficult to obtain con- 
vergence (11,13,14,19). 

Simulations reaching 50 ns also showed ion binding in 
the grooves and found that binding was more extensive for 
potassium than for sodium. It was however noted that dur- 
ing 50 ns, individual ions still only sampled roughly one- 
third of the simulation box, remaining clearly distinguish- 
able from one another, in violation of the assumption of er- 
godicity. In an exceptionally long simulation for the time, 
Perez et al. studied the so-called Dickerson-Drew dode- 
camer (dCGCGAATTCGCG) for 1 |xs in the presence of 
neutralizing sodium ions. High affinity ion binding in the 
AATT minor groove was observed with residence times up 
to 15 ns (20). This study also found a correlation between 
sodium ion binding and minor groove narrowing, following 
up on an earlier controversy in this area surrounding the in- 
terpretation of crystallographic results (1,4). 

When analyzing molecular simulation data, in addition 
to problems linked with convergence, there are questions as 
to how to best measure ion distributions. Deoxyribonucleic 
acid (DNA) is a flexible macromolecule that, given its heli- 
cal conformation, deforms easily by bending, twisting and 
stretching on the timescale of MD simulations. Thus, while 
it is possible to generate an overview of ion distributions by 
superposing the instantaneous DNA conformations from 
successive 'snapshots' (typically spaced in time by 1 ps inter- 
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Figure 1. Left: Schematic view of the curvilinear helicoidal coordinates (CHC). An ion (red dot) is described by a distance D along the curved helical axis 
(black line), a radial distance R from the axis and an angle A from a reference vector which tracks the helical twist of the nucleic acid. At the base pair 
levels, this vector corresponds approximately to the long axis of the base pairs and points toward the 5 f -3' strand. Consequently, A % 90° places an ion 
in the minor groove and A % 270° places an ion in the major groove. Right: Isodensity surfaces (red) for the phosphorus atoms of the AGCT oligomer 
analyzed in the CHC system, then mapped into Cartesian space using the average helical axis of the oligomer (black line). The nucleotides are colored to 
indicate the base sequence (G: blue, C: green, A: red, T: orange). All isodensity plots were obtained using Chimera (39,40). 



vals), the overall motions of the double helix will necessarily 
lead to some 'blurring' of the distributions associated with 
ions that are strongly interacting with DNA. Alternatively, 
one can count ion contacts with specific atoms of DNA us- 
ing a chosen distance cutoff. This gives good local informa- 
tion, but is not well adapted to ions positioned, or fluctuat- 
ing, between clusters of atoms, and it is clearly not adapted 
to analyzing overall distributions. 

We now propose a new approach to analyze the distri- 
bution of ions, or molecules, around DNA based on the 
use of a natural coordinate system for DNA, namely curvi- 
linear helicoidal coordinates (CHC). Since our algorithm, 
Curves +, provides a method for obtaining the helical axis 
of any DNA conformation, it is relatively easy to extend 



this analysis to determine the positions of ions or molecules 
around DNA in terms of their location with respect to the 
helical axis (using the base pair positions to define steps 
along the molecule and the local position of the helical 
grooves). Making this choice implies that overall helical de- 
formations including twisting, bending and stretching can 
be eliminated from the analysis of ion densities. It is then 
possible to plot ion distributions radially from the curvilin- 
ear helical axis (i.e. as a cylindrical distribution function), 
along the length of the helical axis, or around the (curvi- 
linear) helical axis, either for the entire molecule or for any 
chosen zone. We can use this approach to calculate average 
ion populations or local concentrations close to any chosen 
part of the molecule and also analyze the convergence of the 
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Figure 2. Root mean square fluctuations (A) for the phosphorus atoms of the AGCT oligomer are smaller when calculated using the CHC analysis mapped 
into Cartesian space (black lines), than when simply superposing snapshots from the molecular dynamics simulation (red lines). Note that phosphates 
belonging to the two strands are shown consecutively in the 5'-3' direction. 



ion distributions as a function of the duration of the MD 
simulations. Lastly, by using a single, average DNA con- 
formation, we can also map the helicoidal ion coordinates 
back into Cartesian space and obtain 3D density plots that 
do not suffer from the DNA conformational fluctuations in 
the same way as the simple Cartesian superposition of MD 
snapshots and, as discussed below, lead to strikingly clear 
ion distributions. 

This paper describes how curvilinear helicoidal ion pa- 
rameters are defined and calculated. We then apply the new 
method to analyze potassium ion (K + ) and chlorine ion 
(Cl~) distributions around two different DNA oligomers 
that were simulated for 1 (xs in aqueous solution with a 
physiological concentration of KC1. Although K + ions are 
expected to bind more strongly to DNA than Na + be- 
cause of their weaker first hydration shell (15-17,21), the 
results presented here show sequence-specific ion binding 
sites with surprisingly high occupancies. The new method 
enables spatially-localized, time-averaged populations and 
the corresponding molar concentrations to be calculated. 
We show that 1, 2 and 3D graphs using the curvilinear he- 
licoidal coordinates make it easy to understand the overall 
ion distributions. 

The software used for this analysis, 'Canion', is freely 
available, along with the latest version of Curves + and the 
necessary user guides, at the following site: http://bisi.ibcp. 
fr/tools/curves_plus/. We note that Canion can be used not 
only for analyzing data from MD simulations, but can also 
directly read density matrices derived from continuum ap- 
proaches to modeling ion distributions (see, e.g. (22,23)). 
This makes it easy to compare such results with MD simu- 
lations using the curvilinear helicoidal coordinate system. 



MATERIALS AND METHODS 

DNA trajectories 

We illustrate the ion analysis developed here using microsec- 
ond long molecular dynamics trajectories for two 18-mer 
oligomers belonging to the so-called ABC (Ascona B-DNA 
Consortium) database (24-26). Both oligomers are based 
on tetranucleotide repeats, respectively, AGCT and ATCG, 
flanked by GC ends to reduce fraying during the simula- 
tions. Their sequences are GCCTAGCT AGCTA GCTGC 
and GCGCATGC ATGCA TGCGC (where the tetranu- 
cleotide repeat used to refer to the oligomers has been un- 
derlined). 

Each oligomer was simulated using the AMBER pro- 
gram suite (27) and with the current ABC protocol (24). 
This involves the parm99 force field (28,29) with the bscO 
modifications (30). The oligomers were simulated using pe- 
riodic boundary conditions with a truncated octahedral 
cell and a solvent environment consisting of SPC/E water 
molecules (31) (with a minimum thickness of 10 A around 
the solute), leading to 11 621 and 11 546 water molecules 
for AGCT and ATCG, respectively). The DNA net charge 
is neutralized with K + ions and then KC1 ion pairs are added 
to reach a concentration of 1 50 mM KC1. The ions were rep- 
resented using Dang parameters (32). Counterions were ini- 
tially placed at random, at least 5 A from the solute and 3.5 
A from one another. Electrostatic interactions were treated 
with the particle mesh Ewald method (33), using a real- 
space cutoff of 9 A and cubic B-spline interpolation onto a 
charge grid with a 1-A spacing. Lennard- Jones interactions 
were truncated at 9 A. The pair list was built with a buffer re- 
gion and updated whenever a particle moved by more than 
0.5 A. 

The oligomers were initially built with a canonical B- 
DNA conformation. Initial equilibration involved energy 
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Figure 3. Average phosphorus distributions calculated using CHC for the 1 juls AGCT trajectory and plotted in various planes: DA (top left), DR (top 
right), RA (bottom left). The bottom right plot shows the RA plane for the sugar CV atoms from the same trajectory. The blue to red color scale represents 
increasing molarity. DR and RA plots show the average radius of the phosphorus atoms from the helical axis as a white line and as a white circle, respectively. 
The DA and RA plots show the minor groove limits, defined by the average CI' positions, as a white line and as white radial vectors, respectively. RA plots 
also have a vertical radial vector indicating the center of the major groove 

minimization of the solvent, followed by a slow thermal- 
ization. Production simulations were carried out using a 2 
fs time step in an NPT ensemble. The length of chemical 
bonds involving hydrogen were restrained using SHAKE 
(34) and the Berendsen algorithm (35) was used to control 
the temperature and the pressure, with a coupling constant 
of 5 ps. Center of mass motion was removed every ps to limit 
the translational kinetic energy of the solute (36). 

Both oligomers were simulated for a total of 1 (jls and 
conformational snapshots were saved for analysis every 1 
ps (leading to 10 6 data sets per oligomer). 



DNA conformational analysis 

The DNA conformation and the ion distribution were an- 
alyzed directly from the AMBER trajectory output using 
Curves+ (37,38). In the case of DNA, this analysis yields a 
curvilinear helical axis, and a full set of helical, backbone 
and groove parameters that can subsequently be used to 
generate time series or time-averaged probability distribu- 
tions using the program Canal. As discussed below, a new 
option 'ions' is now available in Curves + that triggers the 
calculation of the position of ions, water molecules, or any 
desired solute atoms in the trajectory. This output can then 
be analyzed with the Canion utility program, which is freely 
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Figure 4. Time-averaged K + populations within the DNA grooves for the unique base pair steps (T8pA9, A9pG10, GlOpCll) belonging to the central 
tetranucleotide of the AGCT oligomer for increasing durations (ns) of the molecular dynamics trajectory. 




Figure 5. Averaged conformation of the AGCT oligomer (black line drawing) and the associated average locations of the K + (blue spheres) and CI (green 
spheres) ions calculated for increasing durations of the molecular dynamics trajectory. 
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Figure 6. ID K + distributions averaged over the 1 juls AGCT trajectory. The vertical line in the R plot indicates the radial position of the phosphorus 
atoms, while the shorter distance between the two vertical lines in the A plot delimits the minor groove using the angular position of the sugar CI' atoms 
(see Figure 3). 



available along with Curves +, Canal and other utility soft- 
ware (http://bisi.ibcp.fr/tools/curves_plus/). 



Ion analysis technique 

In order to describe the position of atoms or ions around 
a nucleic acid substrate (or belonging to the substrate it- 
self), we have chosen to use coordinates based on the heli- 
cal axis of the nucleic acid. The Supplementary Information 
to this article contains a detailed description of this proce- 
dure. Briefly, the overall helical axis in Curves + is built from 
a set of helical axis systems consisting of vectors d\ , d2 and 
centered at point y. is the local helical axis and d\di 
are coupled to the helical twist of the molecule, d\ , being 
linked to the long axis of the base pairs. For each ion (or 
atom), we search for the helical axis system where the ion 
lies in the d\d^ plane. If several such systems are found, we 
choose the one leading to the shortest distance between the 
helical axis reference point y and the ion. This distance de- 
fines the radial coordinate R from the ion to the helical axis 



(see Figure 1). The distance of the point 7 along the helical 
axis in units of base pair steps defines the coordinate D (i.e. 
D can vary continuously from 1 to TV within an N bp seg- 
ment) and, finally, the angle of the vector from y to the ion 
defines an angle with respect to a vector d\ which consti- 
tutes the angular coordinate A. Given the definition of d\, 
A ~ 90° corresponds to the center of the minor groove and 
A ~ 270° corresponds to the center of the major groove for 
a canonical B-DNA. 

For a molecular dynamics trajectory (or an ensemble 
of experimental structures, such as that resulting from an 
NMR study), the ion analysis is performed on each snap- 
shot (or experimental structure) using the corresponding 
helical axis and stored in a file. This file is then read by 
the Canion program and the ions positions are accumulated 
in a 3D histogram with a bin size in curvilinear helicoidal 
space of 0.5 A in R, 1/6 of a base pair step (^0.5 A) in D 
and 5° in A. We then analyze the ion distributions in ID (R, 
D or A), in 2D (RA, DA, DR) or in 3D Cartesian space (see 
below). Note that in the case of a 2D RA analysis, we use 
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polar coordinate plots to make the results easier to under- 
stand. Ion distributions can be obtained for the entire space 
surrounding the oligomer, or for any selected zone, defined 
by fixing lower and upper limits on R, D or A. Note that, by 
default, we use an upper limit of R = 30 A since beyond this 
point the solute molecule has little impact on the ion distri- 
bution and the helicoidal coordinate analysis ceases to be 
of interest. 

For any chosen spatial region, we can obtain the time- 
averaged ion populations. However, it is also useful to be 
able to calculate ion densities, or more specifically molari- 
ties. (Note densities in ions. A -3 can be converted to mo- 
larities by dividing by 7V A /10 27 = 6.022 x 10~ 4 , where 10 27 
is the number of cubic angstroms per liter and Na is Avo- 
gadro's number.) Since the volumes of the histogram bins 
in CHC change as a function of R and, if the helical axis 
is curved, also as a function of A at each base pair step, it 
is necessary to calculate the Cartesian volume for each ele- 
ment of the CHC histogram. Note also that curved helical 
axes impose limits on the radial distance at which curvilin- 
ear helicoidal coordinates can be calculated. See the Sup- 
plementary Information for details of these calculations. 

In order to recreate a 3D Cartesian space distribution, the 
helicoidal ion coordinates can be mapped back to Carte- 
sian space with respect to a chosen helical axis, with typi- 
cal choices being the axis of the average nucleic acid struc- 
ture from the corresponding molecular dynamics trajectory, 
or a linear axis corresponding to a regular helix with user- 
defined twist and rise. 

Lastly, when calculating ion densities close to DNA, it 
seems reasonable to exclude the volume occupied by the so- 
lute molecule itself. This has been made possible in the Can- 
ion program by determining which CHC volume elements 
lie inside the van der Waals envelope of the solute molecule 
(using standard Pauling values, P: 1.9 A, C: 1.6 A, N: 1.5 
A, O: 1.4 A, H: 1.2 A). While taking the solute excluded 
volume into account is optional in the Canion program, it 
is particularly useful when comparing the ion molarities in 
the two grooves of helical nucleic acids, and this option has 
been selected in all the results presented here. 

RESULTS AND DISCUSSION 

Before discussing the ion distributions around the AGCT 
and ATGC oligomers, we can first ask whether using the 
CHC system is in fact helpful. To test this we will use the 
phosphorus atoms belonging to the AGCT oligomer as an 
example. 

We start by recording the position of each of these atoms 
from every 1 ps snapshot along the 1 |xs trajectory, using the 
curvilinear helicoidal coordinates defined in the methodol- 
ogy section. Each phosphorus atom is described by its radial 
R, longitudinal D and angular A coordinates (see Figure 1). 
This implies that the most significant overall fluctuations 
of the oligomer (stretching, twisting and bending) will not 
affect the helicoidal coordinates describing its constituent 
atoms. To test this, we first generate an average DNA struc- 
ture from the trajectory. We then calculate its helical axis 
with Curves + (37) and use it to convert the instantaneous 
helicoidal coordinates of the phosphorus atoms from each 



snapshot to a common Cartesian coordinate system. The 
scatter in the positions of each phosphorus atom can then 
be used to generate the isodensity surfaces shown on the 
right of Figure 1 (plotted using Chimera (39,40)). More 
quantitatively, we can calculate root mean square fluctu- 
ations (RMSF) for each phosphorus atom and compare 
the values with those obtained by simply superposing the 
Cartesian coordinates from each snapshot. The results in 
Figure 2 show that the CHC method leads to smaller RMSF 
values in almost all cases. This reduction is particularly in- 
teresting in the center of each strand, where the overall fluc- 
tuations of the oligomer might not have been expected to 
have a significant effect. 

Analyzing the position of DNA atoms in helicoidal co- 
ordinates is also useful for providing reference positions 
to help in analyzing ion distributions. This is illustrated 
in Figure 3 where we plot the phosphorus distributions 
in helicoidal coordinate space using the three possible 2D 
plots: DA (top left), DR (top right) and RA (bottom left). 
The DA plot (obtained by summing the phosphorus den- 
sity over all radial distances 7?) shows an 'unwrapped' and 
'unwound' view of the oligomer. The phosphate densities 
of each strand consequently lead to a column of distribu- 
tions at constant A values, with the smaller separation be- 
tween the two columns corresponding to the minor groove 
of the double helix. In the DR plot (obtained by summing 
over all angles A), all the phosphorus atoms are seen to lie 
at constant distance (R = 10.25 A) from the helical axis. 
This distance can be used to define a useful radius delimit- 
ing the helical grooves, as shown by the white circle in the 
RA plot (obtained by summing over all axial distances D). 
Here, all the phosphorus atom distributions are superposed 
and closely follow the R = 10.25 A radius circle. 

Figure 3 contains one further plot (bottom right) corre- 
sponding to an RA analysis of the sugar CY atoms. These 
atoms constitute a useful guide to the position of the mi- 
nor groove since they reflect the position of the glycosidic 
bonds. Note that the CY densities are tighter than those of 
the phosphorus atoms, since they are less affected by back- 
bone fluctuations, but also because, being closer to the he- 
lical axis, they are affected less by effectively unwinding the 
B-DNA helix in the helicoidal coordinate space (an effect 
more visible for the more distant phosphorus atoms in the 
DA and RA plots). Drawing radial vectors through the cen- 
ter of the CY distributions enables us to define the minor 
groove as lying between A values of 33° and 147°. We take 
remaining values of A to belong to the major groove, since 
we do not attribute regions specifically to the phosphodi- 
ester backbones. We complete the visual reference system 
in RA plots by adding a line along the external bisector of 
the CY vectors to indicate the center of the major groove. 
Note that the CI' vectors can also be used to show the minor 
groove position in ID A plots (where they become vertical 
lines) and, similarly, the phosphorus radius can be added as 
a vertical line in ID R plots (see below). 

We should remark that the usefulness of CHC approach 
depends on having a well-defined helical axis. Curves+ cal- 
culates a curvilinear axis that runs between the terminal 
base pairs of the oligomer. However, if significant end- 
fraying occurs, this can lead to perturbations of the axis 
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Figure 7. 2D K + distributions averaged over the 1 juls AGCT trajectory: DA plane (top left), DR plane (top right), RA plane (bottom). The results are 
plotted as molarities as shown by the color bars, with a blue to red scale indicating increasing values. 



which will then be reflected in the helicoidal coordinates 
used to describe atomic (or ionic) positions. In addition, the 
position of ions beyond the ends of the helical axis will not 
be recorded. 

Before proceeding with the analysis of ion distributions 
around DNA, we need to consider the question of conver- 
gence. As shown below, the main ion densities occur in the 
DNA grooves (R < 10.25 A) and between successive base 
pairs (defined here as TV — 0.2 < D < N + 1.2 for a base 
pair step N V N +1). We have therefore tested convergence 
by looking at the average potassium ion populations as a 
function of time for unique base pair steps belonging to the 
central tetranucleotide and in the associated major or minor 
grooves. The results for the AGCT oligomer are shown in 
Figure 4 and the corresponding data for ATGC in Supple- 
mentary Figure S4. Both these figures show that at least 300 
ns of simulation are necessary to achieve stable ion popula- 
tions (and, in some cases, non-negligible changes can occur 
up to 500 ns). 



A second test of convergence is shown in Figure 5, where 
we have simply generated average Cartesian coordinates for 
the AGCT oligomer using 1 ps snapshots and including the 
position of all the K + (blue) and Cl~ (green) ions. When the 
average is made over a few tens of nanoseconds of simula- 
tion, the ions tend to sample limited volumes of the simula- 
tion cell (15) and are thus clearly distinguishable. However, 
after 1 |xs of simulation the average ion positions all coin- 
cide at the origin of the simulation cell, showing they have 
had time to thoroughly sample the full simulation cell in an 
equivalent manner, satisfying the principle of ergodicity. 

Having shown that reasonable convergence has been 
achieved, we turn to an overall analysis of the potassium ion 
distributions around the AGCT oligomer. We begin with 
ID R, D and A plots shown in Figure 6. The radial distance 
R plot, that is also a radial distribution function, shows 
two peaks within the grooves and one beyond the phos- 
phorus atom radius (indicated by the vertical line). Note 
that the maximum molarity for the peak closest to the he- 
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Figure 8. 3D K + distributions obtained by mapping CHC analysis of the 1 fxs AGCT trajectory into Cartesian space with respect to the average DNA 
structure, shown as a line drawing on the left (G: blue, C: green, A: red, T: orange) and as a gray solvent accessible surface on the right. Molarity isodensity 
surfaces for potassium ions are plotted at 15 M (solid red) and 5 M (green mesh). 



lical axis reaches almost 14 times the bulk value, whereas, 
at 30 A the value has fallen slightly below the effective bulk 
molarity (0.336 M, taking into account the extra K + ions 
added to achieve neutrality). The axis distance D plot shows 
two interesting features: first, three repeated profiles that 
are in register with the repeating base sequence and which 
show strong molarity peaks at three consecutive GpC steps; 
second, a decrease in ion density toward the ends of the 
oligomer, reflecting the decreasing phosphate density (41). 
Lastly, the angular A plot shows that potassium ion density 
peaks occur in both grooves. 

We can further refine this analysis by passing to 2D repre- 
sentations. Figure 7 shows the DA and RA plots for potas- 
sium ions, again around the AGCT oligomer (with radial 
sampling limited to 15 A). This enables us to identify the R 
peak closest to the helical axis (seen in Figure 6) as ion bind- 
ing in the center of the major groove, while on the minor 
groove side there is a more diffuse density with two max- 



ima (see the DA and RA plots). Before filling in the details 
of the potassium ion distributions, it is worth contrasting 
these results with those for the chlorine anions. Supplemen- 
tary Figure S5 shows that chlorine anions have almost no 
density within the phosphorus radius, although they begin 
to penetrate when the phosphate density decreases toward 
the ends of the oligomer (see the DR plot). As expected, this 
penetration mainly occurs on the side of the wider major 
groove (see the RA plot). 

In order to complete the potassium ion analysis, it is use- 
ful to look at the 3D ion distributions reconstructed from 
the instantaneous helicoidal coordinates using a single, av- 
erage helical axis (as for the RMSF calculations in Figure 
2). These distributions are shown as isodensity plots in Fig- 
ure 8. The left image shows the average DNA conformation 
as a wire diagram with K + ion distributions at two isoden- 
sity levels: 15 M (red solid volumes) and 5 M (green mesh). 
The right image shows the same isodensities, but with a sur- 
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Figure 9. 2D RA K + distributions for various dinucleotide steps within the central tetranucleotides of the AGCT and ATGC oligomers, averaged over 
the corresponding 1 juus trajectories. The blue to red color scale indicates increasing molarity values. In each plot, the numbers in the upper and lower 
semicircles indicate the average K + occupancy in the major and minor grooves, respectively. 
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Figure 10. Variations in K + molarity along the DNA grooves (solid lines) of the AGCT oligomer. Values are averaged over the 1 juls AGCT trajectory and 
compared with variations in groove width (dotted lines). Groove width variations are plotted in A with respect to the respective minimal values (major: 
10.0 A, minor: 3.5 A) on the same scale as the molarities. Left: major groove. Right: minor groove. 
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face view of DNA that helps to locate the ions within the 
major and minor grooves. As already seen in Figure 4, the 
largest ion densities occur in the major groove of the GpC 
steps, whereas the largest minor groove densities occur at 
TpA steps. 

At this stage, it is useful to compare the AGCT oligomer 
with the results obtained for the ATGC sequence. Although 
the latter oligomer has a similar overall B-DNA structure 
and the same 50/50 balance of AT and GC base pairs, its 
ion distribution is rather different as shown by the 3D im- 
ages in Supplementary Figure S6. Using the same isodensity 
surfaces as Figure 8 makes it clear that local K + densities oc- 
cupy smaller volumes in both grooves. In the major groove, 
they are again localized at GpC, while in the minor groove 
the densities are considerably weaker and associated with 
ApT steps. 

Using the 3D information we can now return to 2D RA 
plots for a detailed and quantitative analysis for individ- 
ual base pair steps and for each groove. Figure 9 shows 
the results for unique base pair steps of the central tetranu- 
cleotide of each oligomer: TpA, ApG and GpC for the 
AGCT oligomer and ApT, TpG and GpC for the ATGC 
oligomer. In each case, the average ion populations within 
the grooves are shown (corresponding to the data already 
plotted in Figure 4). It is striking that GpC steps in the 
AGCT oligomer have potassium ions strongly localized in 
their major grooves almost 80% of the time. This value de- 
creases somewhat, but still exceeds 60%, for the same step 
in the ATGC oligomer. These main ion-binding sites lead to 
local molarities above 50 M, almost 150 times higher than 
the bulk potassium ion concentrations. Secondary ion bind- 
ing sites, occupied roughly 25-35% of the time, occur in the 
minor groove of TpA steps and the major groove of ApG 
steps for AGCT, and also in the major grooves of ApT and 
TpG steps for ATGC. 

Finally, we consider the question of whether cation bind- 
ing within the grooves of DNA is related to fluctuations in 
groove width. This question has been the subject of exten- 
sive discussions (1,4,12,15,20), it being argued that bound 
cations should be able to counteract phosphate repulsion 
and lead to groove narrowing. Conclusions were however 
unclear, in part because it was difficult to define what 'be- 
ing within the grooves' meant from a geometrical point of 
view. We can now revisit this question in the light of the 
CHC ion analysis. For this, we have added the groove width 
variations to the ID axis distance D plots of K + ion den- 
sity, where the summation over R has now been limited to 
the grooves (R = 0-10.25 A) and the summation of the an- 
gles A has been limited to either the minor or major groove. 
Figure 10 shows the results for AGCT oligomer. For this 
case, there is a clear oscillation of both groove widths (dot- 
ted lines), in register with the tetranucleotide repeating se- 
quence. This oscillation is related to peaks in the ion density, 
but note that the grooves are wide when the ion density in 
the grooves is high. These findings are contradictory to an 
argument based on groove narrowing due to the compen- 
sation of phosphate repulsion. In the ATGC oligomer (see 
Supplementary Figure S7), the groove widths vary signifi- 
cantly less along the sequence (other than at the ends of the 
oligomer), but again there is no indication that ion bind- 



ing leads to groove narrowing. Confirming these findings 
will however require the study of a wider range of base se- 
quences. 

CONCLUSION 

We propose a novel method for analyzing the distributions 
of ions or molecules (notably solvent molecules) around he- 
lical nucleic acids. By using curvilinear helicoidal coordi- 
nates to locate the particles we eliminate the 'blurring' of 
the distributions for ions whose movements are coupled to 
the large amplitude fluctuations of the double helix, notably 
bending, stretching and twisting. In addition, the helicoidal 
coordinates make it possible to calculate ion populations, or 
densities, in any zone surrounding the solute molecule. One 
can visualize the resulting ion distributions in helicoidal co- 
ordinates in ID or 2D and also map the helicoidal coordi- 
nates back into Cartesian space to visualize 3D Cartesian 
space ion densities. Applying this technique to the analy- 
sis of 1 fjus molecular dynamics trajectories of two DNA 
oligomers enables us to quantify the convergence of the 
ion distributions (which requires at least 300 ns of simula- 
tion) and to identify strongly localized, sequence-dependent 
ion binding sites (without requiring the use of distance cut- 
offs with respect to specific solute atoms). This analysis can 
be applied to any helical nucleic acid conformations (de- 
rived from molecular simulations or from experiment) and 
to any surrounding species such as ions, water molecules 
or ligands. The software necessary to carry out the analy- 
sis (Curves + and associated utility programs) is freely avail- 
able. 
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